procedure doall()
#(ssn_prefix, ssn_value,
#                aamin, aamax, aastep, 
#                ffmin, ffmax, ffstep)

char    ssn_prefix   {"SN",prompt="Base name for directories "}
int     ssn_value    {100,min=1,max=1000,prompt="SN"}
int     aamin        {1,min=1,max=10,prompt="Minimum a"}
int     aamax        {5,min=1,max=10,prompt="Maximum a"}
int	aastep       {1,min=1,max=10,prompt="a step"}
int     ffmin        {1,min=1,max=10,prompt="Minimum fwhm"}
int     ffmax        {5,min=1,max=10,prompt="Maximum fwhm"}
int	ffstep       {1,min=1,max=10,prompt="fwhm step"}
int	start	     {1,prompt="first cont map"}

begin

char sn_prefix       # prefix for directory names
int sn_value
int amin,amax,astep  # spacing variables
int fmin,fmax,fstep  # fwhm variables
int mc               # counts MCarlo loop

int a, f     # spacing and fwhm - loops
char new_dir
int nmask    # reads file database/Nmask to know how many apertures were found
int na, nf   # number of points in spacing and fwhm
char dimstr  # to be used by wtext when building the contamination image
char title   # used by contour to create a title for the plot

sn_prefix = ssn_prefix
sn_value = ssn_value

amin = aamin
amax = aamax
astep = aastep
fmin = ffmin
fmax = ffmax
fstep = ffstep

# should check, this should be an integer
# a - spacing - number of different spacings
na = (amax-amin)/astep
# f - fwhm - number of different fwhms
nf = (fmax-fmin)/fstep
# contamination, quick look
  if ( access("cont.ql") ) {
    delete("cont.ql",verify=no)
  }
  !touch cont.ql
# contamination, linear gauss fit   print ("LA VAI")
  if ( access("cont.lin") ) {
    delete("cont.lin",verify=no)
  }
  !touch cont.lin

for ( mc = start ; mc<start+1 ; mc=mc+1 ) { # Monte Carlo index
  for ( a=amin ; a<=amax ; a=a+astep ){ # loops in fiber spacing
    for ( f=fmin ; f<=fmax ; f=f+fstep ){ # loops in fwhm
# the aperture width is the same as the fiber spacing
      print("Looping fibers")
	apresize.llimit = -1 * a / 2.0 
	apresize.ulimit =  a / 2.0 
	print("Limits = ",apresize.llimit,apresize.ulimit)
	new_dir = sn_prefix // str(sn_value) // "_a" // str(a) // ".0" // "_f" // str(f) // ".0"
	chdir (new_dir)
	print (new_dir)
	if (imaccess("G8IV"))
	    imdelete("G8IV")
	imcopy("G8IV."//mc,"G8IV",verbose=no)
	!cleansimul
	mask
	nmask = 10
	scan(nmask,<"database/Nmask")
# this nmask==25 test
# is used to discard bad fits.  if I could not find all 25 fibers 
# then either SN is too small or the fibers are too packed.
	if (nmask == 25 )	{ 
# print("CALFIB")
	  calfib(nome="mask")
# here I should loop around all object files apertures first
# print("QL")
	  ifu (image="G8IV",ql=yes)
# here I should loop around all object files then linear gauss fit
	  print("LGAUSS")
	  ifu (image="G8IV",ql=no)
	  wtext (input="G8IV.lin.ldisp.fits",output="G8IV.lin.ldisp.txt")
	  wtext (input="G8IV.ql.ldisp.fits",output="G8IV.ql.ldisp.txt")
	  !cont < G8IV.lin.ldisp.txt >> ../cont.lin
	  !cont < G8IV.ql.ldisp.txt >> ../cont.ql
	}
	else {
	  print "ECHO 1.00"
	    !echo "1.0000" >> ../cont.lin
	    !echo "1.0000" >> ../cont.ql
	    }
      chdir ("..")
	}
  }
  
  print ("Done calculating, creating contamination maps")
    dimstr = str(na+1) // "," // str(nf+1)
    
    print("Dimstr = ",dimstr)
    
  print("RTEXT1")
  rtext(input="cont.lin",output="cont.lin",otype="",header=no,pixels=yes,
	nskip=0,dim=dimstr)
  print("RTEXT2")
  rtext(input="cont.ql",output="cont.ql",otype="",header=no,pixels=yes,	nskip=0,dim=dimstr)
    
#title = "Aperture extraction - y spacing - x fibers width in sigma" \
#        // sn_prefix // str(sn_value)
    contour(image="cont.lin",
	    floor=0.1,ceiling=0.8,zero=0.,
	    ncontours=4,interval=0.2,nhi=-1,dashpat=528,device="stdgraph",
	    title="A",
	    preserve=yes,label=yes,fill=no,xres=64,yres=64,perimeter=yes,
	    vx1=0.,vx2=0.,vy1=0.,vy2=0.,subsample=no,append=no,mode="ql")
    
#title = "Gaussian extraction - y spacing - x fibers width in sigma" \
#        // sn_prefix // str(sn_value)
    contour(image="cont.ql",floor=0.1,ceiling=0.8,zero=0.,
	    ncontours=4,interval=0.2,nhi=-1,dashpat=528,device="stdgraph",
	    title="B",
	    preserve=yes,label=yes,fill=no,xres=64,yres=64,perimeter=yes,
	    vx1=0.,vx2=0.,vy1=0.,vy2=0.,subsample=no,append=no,mode="ql")
    print("RENAME")
    rename("cont.lin","cont.lin."//mc)
    rename("cont.ql","cont.ql."//mc)
    !touch cont.lin cont.ql
    imrename("cont.lin","cont.lin."//mc)
    imrename("cont.ql","cont.ql."//mc)
 }
    end
    
